#include "base.h"

int main(int argc,char **argv)
{

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@///

	  //////////////////
	  // Input   ///////
	  //////////////////
          char kappa_sea[1024];
          int conf_ini, jump;
          int  Nex;
          int Nmass,Nr;
          int NOa2;
	  int Nf;
          double beta, Uplaq;
          char dir[1024];
          char gauge_fixing[1024];
          char gauge_action[1024];



          ////// Parsing command line

          for(int i=1; i<argc; i++) {
            if(strncmp(argv[i],"-kappa_sea",10)==0)  sscanf(argv[i],"-kappa_sea=%s",kappa_sea);
            if(strncmp(argv[i],"-conf_ini",9)==0)  sscanf(argv[i],"-conf_ini=%d",&conf_ini);
            if(strncmp(argv[i],"-jump",5)==0)  sscanf(argv[i],"-jump=%d",&jump);
            if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
            if(strncmp(argv[i],"-Nex",4)==0)  sscanf(argv[i],"-Nex=%d",&Nex);
            if(strncmp(argv[i],"-Nmass",6)==0)  sscanf(argv[i],"-Nmass=%d",&Nmass);
            if(strncmp(argv[i],"-Nr",3)==0)  sscanf(argv[i],"-Nr=%d",&Nr);
	    if(strncmp(argv[i],"-Nf",3)==0)  sscanf(argv[i],"-Nf=%d",&Nf);
            if(strncmp(argv[i],"-dir",4)==0)  sscanf(argv[i],"-dir=%s",dir);
            if(strncmp(argv[i],"-gauge_fixing",13)==0)  sscanf(argv[i],"-gauge_fixing=%s",gauge_fixing);
            if(strncmp(argv[i],"-gauge_action",13)==0)  sscanf(argv[i],"-gauge_action=%s",gauge_action);
            if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T);
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);
            if(strncmp(argv[i],"-beta",5)==0)  sscanf(argv[i],"-beta=%lf",&beta);
            if(strncmp(argv[i],"-Uplaq",6)==0)  sscanf(argv[i],"-Uplaq=%lf",&Uplaq);
            if(strncmp(argv[i],"-NOa2",5)==0)  sscanf(argv[i],"-NOa2=%d",&NOa2);

          }

   /////////////////////////
  //// Momenta set
  //////////////////////////////
  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom.dat");
  read_momenta_data(momentafile);

	  //////////////////////////////////
          // IO files
          ///////////////////////////////

	 char file[1024];
	 char z_output[1024];
	 char dir_qav[1024];	 
         sprintf(dir_qav,"%s/QAV",dir);

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

	//////////////////////////////////
	// Read Zs
       // NOTE: Lambda & Zq MUST be reals (check in the last step)
       // I work with reals from now on
	/////////////////////////////

        printf("====================\n");
	printf(" READ Z             \n");
	
        double ****Zq=malloc(sizeof(double***)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust]=malloc(sizeof(double**)*(Nmass));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) Zq[iclust][imass]=malloc(sizeof(double*)*(Nmom));
	for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) Zq[iclust][imass][imom]=malloc(sizeof(double)*(Nr));

        double *****Lambda=malloc(sizeof(double****)*(16));
        for (int igamma=0; igamma<16 ; igamma++) Lambda[igamma]=malloc(sizeof(double***)*(Nclust+1));
        for (int igamma=0; igamma<16 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) Lambda[igamma][iclust]=malloc(sizeof(double**)*(Nmass*Nmass));
        for (int igamma=0; igamma<16 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass*Nmass; imass++) Lambda[igamma][iclust][imass]=malloc(sizeof(double*)*(Nmom));
	for (int igamma=0; igamma<16 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass*Nmass; imass++) for (int imom=0; imom<Nmom; imom++) Lambda[igamma][iclust][imass][imom]=malloc(sizeof(double)*(Nr));


      double **aux_allmomenta=malloc(sizeof(double*)*(Nmom));
     for (int imom=0; imom<Nmom; imom++) aux_allmomenta[imom]=malloc(sizeof(double)*(Nclust+1));

    // Zq
    printf("Zq\n");
    for (int r=0; r<2; r++)
    for (int imass=0; imass<Nmass; imass++)
    {
         sprintf(file,"Zq_r%d_m%d",r,imass);

         read_input_momentabymomenta(dir,file,aux_allmomenta);
	 for (int imom=0; imom<Nmom; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust][imass][imom][r]=aux_allmomenta[imom][iclust];
    }

    // Lambda 
   printf("ZO\n");
   char gammaprint[8];
   for (int igamma=0; igamma<16; igamma++)
   {
        whichGamma(igamma,gammaprint);
	printf("%s\n",gammaprint);
        for (int r=0; r<2; r++)
	for (int im2=0; im2<Nmass; im2++) for (int im1=0; im1<Nmass; im1++)
        {
         int imass=im2*Nmass+im1;
	 sprintf(file,"Lambda%s_r%d_m%d%d",gammaprint,r,im1,im2);
         read_input_momentabymomenta(dir,file,aux_allmomenta);
         for (int imom=0; imom<Nmom; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) Lambda[igamma][iclust][imass][imom][r]=aux_allmomenta[imom][iclust];


       }
    }

	
     printf("====================\n\n\n");

     //////////////////////////////////////////////////////////
     // Promedium r0 r1
     /////////////////////////////////////////////////////////////


     printf("====================\n");
     printf(" PROMEDIUM r0 r1    \n");

     double ***Zqr=malloc(sizeof(double**)*(Nclust+1));
     for (int iclust=0; iclust<Nclust+1; iclust++) Zqr[iclust]=malloc(sizeof(double*)*(Nmass));
     for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) Zqr[iclust][imass]=malloc(sizeof(double)*(Nmom));


      double ****Lambdar=malloc(sizeof(double***)*(16));
      for (int igamma=0; igamma<16 ; igamma++) Lambdar[igamma]=malloc(sizeof(double**)*(Nclust+1));
      for (int igamma=0; igamma<16 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) Lambdar[igamma][iclust]=malloc(sizeof(double*)*(Nmass*Nmass));
      for (int igamma=0; igamma<16 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass*Nmass; imass++) Lambdar[igamma][iclust][imass]=malloc(sizeof(double)*(Nmom));

      for (int iclust=0; iclust<Nclust+1; iclust++) for (int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) 
      		Zqr[iclust][imass][imom]=(Zq[iclust][imass][imom][0]+Zq[iclust][imass][imom][1])/2.;

     for (int igamma=0; igamma<16 ; igamma++)for (int iclust=0; iclust<Nclust+1; iclust++) for (int imass=0; imass<Nmass*Nmass; imass++) for (int imom=0; imom<Nmom; imom++) 
	{
		 Lambdar[igamma][iclust][imass][imom]=(Lambda[igamma][iclust][imass][imom][0]+Lambda[igamma][iclust][imass][imom][1])/2.; 
	}
  


     printf("====================\n\n\n");


     /////////////////////////////////
     // Compute LambdaG where:
     ////// LambdaS=LambdaS0
     /////  LambdaV=LambdaV1+LambdaV2+LambdaV3+LambdaV4
     /////  LambdaP=LambdaP5
     /////  LambdaT=LambdaT1+LambdaT2+LambdaT3+LambdaB1+LambdaB2+LambdaB3
     /////  LambdaA=LambdaA1+LambdaA2+LambdaA3+LambdaA4
     ////////////////////////////////////

     printf("============================\n");
     printf(" COMPUTE ZS SV ZP ZT ZA \n");

      double ****LambdaGr=malloc(sizeof(double***)*(5));
      for (int igamma=0; igamma<5 ; igamma++) LambdaGr[igamma]=malloc(sizeof(double**)*(Nclust+1));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) LambdaGr[igamma][iclust]=malloc(sizeof(double*)*(Nmass*Nmass));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass*Nmass; imass++) LambdaGr[igamma][iclust][imass]=malloc(sizeof(double)*(Nmom));


      for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass*Nmass; imass++) for (int imom=0; imom<Nmom; imom++) 
      {
	// LambdaS
	LambdaGr[0][iclust][imass][imom]=Lambdar[0][iclust][imass][imom];
	// LambdaV
	LambdaGr[1][iclust][imass][imom]=(Lambdar[1][iclust][imass][imom]+Lambdar[2][iclust][imass][imom]+Lambdar[3][iclust][imass][imom]+Lambdar[4][iclust][imass][imom])/4.;
	// LambdaP
	LambdaGr[2][iclust][imass][imom]=Lambdar[5][iclust][imass][imom];
	// LambdaT
	LambdaGr[3][iclust][imass][imom]=(Lambdar[6][iclust][imass][imom]+Lambdar[7][iclust][imass][imom]+Lambdar[8][iclust][imass][imom]+Lambdar[9][iclust][imass][imom]+Lambdar[10][iclust][imass][imom]+Lambdar[11][iclust][imass][imom])/6.;
	// LambdaA
	LambdaGr[4][iclust][imass][imom]=(Lambdar[12][iclust][imass][imom]+Lambdar[13][iclust][imass][imom]+Lambdar[14][iclust][imass][imom]+Lambdar[15][iclust][imass][imom])/4.;

      }	
	
      double ****ZGr=malloc(sizeof(double***)*(5));
      for (int igamma=0; igamma<5 ; igamma++) ZGr[igamma]=malloc(sizeof(double**)*(Nclust+1));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) ZGr[igamma][iclust]=malloc(sizeof(double*)*(Nmass*Nmass));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass*Nmass; imass++) ZGr[igamma][iclust][imass]=malloc(sizeof(double)*(Nmom));


   for (int igamma=0; igamma<5; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for (int imom=0; imom<Nmom; imom++)
    for (int im2=0; im2<Nmass; im2++) for (int im1=0; im1<Nmass; im1++)
    {
        int imass=im2*Nmass+im1;
        ZGr[igamma][iclust][imass][imom]=sqrt(Zqr[iclust][im1][imom]*Zqr[iclust][im2][imom])/LambdaGr[igamma][iclust][imass][imom];
    }


    printf("====================================================\n\n\n");

     //////////////////////////////////////////////
     //     EQUIVALENT MOMENTA
     ///////////////////////////////////////////////
     printf("==================================\n");
     printf(" PROMEDIUM OVER EQUIVALENT SINP2 \n");


     Equivalent_Momenta(SinP2);
        //for (int imom=0; imom<Nmom_eq; imom++) printf("%d %8.5f \n", imom, SinP2_eq[imom]);


      double ****ZGr_eq=malloc(sizeof(double***)*(5));
      for (int igamma=0; igamma<5 ; igamma++) ZGr_eq[igamma]=malloc(sizeof(double**)*(Nclust+1));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) ZGr_eq[igamma][iclust]=malloc(sizeof(double*)*(Nmass*Nmass));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass*Nmass; imass++) ZGr_eq[igamma][iclust][imass]=malloc(sizeof(double)*(Nmom_eq));

     double ***Zqr_eq=malloc(sizeof(double**)*(Nclust+1));
     for (int iclust=0; iclust<Nclust+1; iclust++) Zqr_eq[iclust]=malloc(sizeof(double*)*(Nmass));
     for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) Zqr_eq[iclust][imass]=malloc(sizeof(double)*(Nmom_eq));

     double *aux_mom=malloc(sizeof(double)*Nmom);
     double *aux_eq_mom=malloc(sizeof(double)*Nmom_eq); 

   for (int igamma=0; igamma<5; igamma++)
     for (int im2=0; im2<Nmass; im2++) for (int im1=0; im1<Nmass; im1++)
     {
       int imass=im2*Nmass+im1;	
       for (int iclust=0; iclust<Nclust+1; iclust++)
       {
         for (int imom=0; imom<Nmom; imom++) aux_mom[imom]=ZGr[igamma][iclust][imass][imom] ;
         Average_over_Equivalent_Momenta(SinP2,aux_mom,aux_eq_mom);
         for (int imom=0; imom<Nmom_eq; imom++)  ZGr_eq[igamma][iclust][imass][imom]=aux_eq_mom[imom];
       }
     }
    
   for (int imass=0; imass<Nmass; imass++)     
   for (int iclust=0; iclust<Nclust+1; iclust++)
   {
    for (int imom=0; imom<Nmom; imom++) aux_mom[imom]=Zqr[iclust][imass][imom];
    Average_over_Equivalent_Momenta(SinP2,aux_mom,aux_eq_mom);
    for (int imom=0; imom<Nmom_eq; imom++)  Zqr_eq[iclust][imass][imom]=aux_eq_mom[imom];

   }


     printf("============================\n\n\n");

///***************************************************************************************************************///////
//********************************************************************************************************************///
///***************************************************************************************************************///////

    // write output

      printf(" WRITING OUTPUT IN %s\n",dir);

     double **aux_allmomenta_real=malloc(sizeof(double*)*(Nmom));
     for (int imom=0; imom<Nmom; imom++) aux_allmomenta_real[imom]=malloc(sizeof(double)*(Nclust+1));

     double *aux_clust_real=malloc(sizeof(double)*(Nclust+1));

    ////////////////
    // Zq R-PROM
    ////////////////////
    for (int imass=0; imass<Nmass; imass++)
     {
         sprintf(z_output,"%s/Z_ALLP/Zq_m%d.out",dir,imass);
         for (int imom=0; imom<Nmom; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta_real[imom][iclust]=Zqr[iclust][imass][imom];
         write_output_allmomenta_real(z_output,aux_allmomenta_real,SinP2,Nmom);

         sprintf(file,"Zq_m%d",imass);
         write_output_momentabymomenta_real(dir,file,aux_allmomenta_real,SinP2,Nmom);
     }

	// QAVERAGE
    for (int imass=0; imass<Nmass; imass++)
     {

         sprintf(file,"Zq_m%d",imass);
	 for (int imom=0; imom<Nmom_eq; imom++){
	  for (int iclust=0; iclust<Nclust+1; iclust++) 
		aux_clust_real[iclust]=Zqr_eq[iclust][imass][imom];
	  write_output_onemomenta_real(dir_qav,file,aux_clust_real,SinP2_eq,imom);
    	 } 

         sprintf(file,"%s/Z_ALLP/Zq_m%d.dat",dir_qav,imass);
         FILE *fout=fopen(file,"w");
         if(fout == (FILE *)NULL) {
              fprintf(stderr,"Unable to open file %s for writing\n",file);
              return EXIT_FAILURE;
         }
         for (int imom=0; imom<Nmom_eq; imom++) { 
		for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust_real[iclust]=Zqr_eq[iclust][imass][imom]; 
		fprintf(fout,"%8.10f %8.10f %8.10f  \n",SinP2_eq[imom],Zqr_eq[Nclust][imass][imom],error_jack_real(aux_clust_real,Nclust));
	  }
      fclose(fout);
     }
     /////////////////////
    // Lambda_OPROM R-PROM 
    //////////////////// 


   for (int igamma=0; igamma<5; igamma++)
   {
        whichGammaProm(igamma,gammaprint);
        for (int im2=0; im2<Nmass; im2++) for (int im1=0; im1<Nmass; im1++)
         {
         int imass=im2*Nmass+im1;

         for (int imom=0; imom<Nmom; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta[imom][iclust]=LambdaGr[igamma][iclust][imass][imom];

         sprintf(file,"Lambda%s_m%d%d",gammaprint,im1,im2);
         write_output_momentabymomenta_real(dir,file,aux_allmomenta,SinP2,Nmom);


          }
   }

   ///////////////////////////////
   /// ZO R-PROM & GAMMA-PROM
   ///////////////////////////////////
   for (int igamma=0; igamma<5; igamma++)
   {
        whichGammaProm(igamma,gammaprint);
        for (int im2=0; im2<Nmass; im2++) for (int im1=0; im1<Nmass; im1++)
         {
         int imass=im2*Nmass+im1;

         for (int imom=0; imom<Nmom; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta_real[imom][iclust]=ZGr[igamma][iclust][imass][imom];
         sprintf(file,"Z%s_m%d%d",gammaprint,im1,im2);
         write_output_momentabymomenta_real(dir,file,aux_allmomenta_real,SinP2,Nmom);


          }
   }

   // QAVERAGE
   for (int igamma=0; igamma<5; igamma++)
   {
        whichGammaProm(igamma,gammaprint);
        for (int im2=0; im2<Nmass; im2++) for (int im1=0; im1<Nmass; im1++)
         {
         int imass=im2*Nmass+im1;

	 sprintf(file,"Z%s_m%d%d",gammaprint,im1,im2);
         for (int imom=0; imom<Nmom_eq; imom++) {
		for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust_real[iclust]=ZGr_eq[igamma][iclust][imass][imom];
		 write_output_onemomenta_real(dir_qav,file,aux_clust_real,SinP2_eq,imom);
	}
          sprintf(file,"%s/Z_ALLP/Z%s_m%d%d.dat",dir_qav,gammaprint,im1,im2);
          FILE *fout=fopen(file,"w");
          if(fout == (FILE *)NULL) {
              fprintf(stderr,"Unable to open file %s for writing\n",file);
              return EXIT_FAILURE;
          }
         for (int imom=0; imom<Nmom_eq; imom++) { 
		for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust_real[iclust]=ZGr_eq[igamma][iclust][imass][imom];
		fprintf(fout,"%8.10f %8.10f %8.10f  \n",SinP2_eq[imom],ZGr_eq[igamma][Nclust][imass][imom],error_jack_real(aux_clust_real,Nclust));
	}
      fclose(fout);


          }
   }

   return 0;

}
